Mark D. Scheuerell
Fish Ecology Division, Northwest Fisheries Science Center, National Marine Fisheries Service, National Oceanic and Atmospheric Administration, Seattle, WA USA, mark.scheuerell@noaa.gov

0.0.1 Version

This is version 0.17.08.01.


1 Initialization

library(readr)
library(reshape2)
library(MARSS)

2 Data

For these analyses we will use time series of counts of adult spring/summer Chinook salmon from the Snake River basin in Oregon and Idaho. This group of fish forms an Evolutionarily Significant Unit, which was listed as threatened under the US Endangered Species Act in 1992.

The data were compiled largely through the efforts of Tom Cooney (Northwest Fisheries Science Center, National Marine Fisheries Service, National Oceanic and Atmospheric Administration, Portland, OR USA, tom.cooney@noaa.gov).

2.1 Assets

The assets data file has the following columns of information:

  • pop = abbreviated name of the population
  • code = 5-letter abbreviation for MPG (2 letters) & population (3 letters)
  • mpg = name of 1 of 5 Major Population Groups
  • ha = estimated area (hectares) of available spawning habitat
  • category = hatchery supplementation category (“impact” or “control”)
  • brood_yr = year when spawning occurred
  • nS = total number of spawning fish
  • phos = proportion of hatchery-origin spawners
rawdata <- read.csv("https://raw.githubusercontent.com/mdscheuerell/CAPM/master/data/srss_chin.csv")
head(rawdata)
##      pop  code                   mpg      ha category brood_yr   nS phos
## 1 CathCr GRCAT Grande Ronde / Imnaha 29.1522   impact     1955  122    0
## 2 CathCr GRCAT Grande Ronde / Imnaha 29.1522   impact     1956 2606    0
## 3 CathCr GRCAT Grande Ronde / Imnaha 29.1522   impact     1957 1754    0
## 4 CathCr GRCAT Grande Ronde / Imnaha 29.1522   impact     1958  486    0
## 5 CathCr GRCAT Grande Ronde / Imnaha 29.1522   impact     1959  712    0
## 6 CathCr GRCAT Grande Ronde / Imnaha 29.1522   impact     1960 3161    0

The first thing to do is remove hatchery-origin spawners from the counts, and then convert the counts to log-density based on the estimated spawning area for each population. In a couple of years, incomplete censuses with small sample sizes led to estimated fractions of wild fish equal to zero. I will treat those as NA.

dat <- rawdata[,c("mpg", "code", "pop", "category", "brood_yr")]
## calculate log-density
dat$Sw <- round(log(rawdata$nS*(1-rawdata$phos)/rawdata$ha), 2)
## replace possible -Inf with NA
dat$Sw[is.infinite(dat$Sw)] <- NA

For easier bookkeeping, I’ll also switch the Imnaha’s MPG over to the Grande Ronde.

## change Imnaha name for easier sorting later
dat[,"code"] <- sub("IRMAI", "GR_IR", dat[,"code"], fixed=TRUE)
## MPG abbreviations
mpg_ID <- c("GR", "MF", "SF", "SR")

Only a few of the populations have observations going all the way back to the 1950s, so I’ll trim the data accordingly.

## use generation time of 4 years
gt <- 4
## first brood year to consider
yr_first <- 1960
## last brood year to consider
yr_last <- 2011
## years of interest
t_index <- seq(yr_first, yr_last-gt)
## length of time period
TT <- length(t_index)
## subset data
dat <- subset(dat, brood_yr>=yr_first & brood_yr<=yr_last)

I will use the MARSS package to fit the CAPM model, which requires the data to be formatted as an [N x T] matrix.

## reshape assets
assets <- dcast(dat, brood_yr ~ code, value.var="Sw")[,-1]
## compute differences
assets <- t(apply(assets, 2, diff, lag=gt))
matplot(t(assets), type="b")

## number of popns
NN <- dim(assets)[1]
## names of popns
names_pop <- rownames(assets)

Finally, here are some nicer, longer names for the various populations.

# set up naming & ordering scheme
longnames <- c(
  "Wenaha",
  "Grande Ronde",
  "Catherine",
  "Minam",
  "Lostine",
  "Imnaha",
  "SF Salmon",
  "Secesh",
  "EFSF Salmon",
  "Big",
  "Sulphur",
  "Bear Valley",
  "Marsh",
  "Loon",
  "Camas",
  "Yankee Fork",
  "Valley",
  "UM Salmon",
  "LM Salmon",
  "EF Salmon",
  "Lemhi")
name_ord <- c(6,3,5,4,2,1,12,10,15,14,13,11,9,7,8,20,21,19,18,17,16)
names_mat <- data.frame(ord=name_ord, ID=names_pop)
names_mat <- names_mat[order(names_mat$ord),]
names_mat$long <- longnames

2.2 Market index

The first “market”" index is based on the monthly North Pacific Gyre Oscillation (NPGO) data provided by Emanuele Di Lorenzo, which are available here. We begin by downloading the raw NPGO data and viewing the metadata.

## URL for NPGO data
url_NPGO <- "http://www.o3d.org/npgo/npgo.php"
## raw NPGO data 
all_NPGO <- read_lines(url_NPGO)
## line with data headers
hdr_NPGO <- which(lapply(all_NPGO,grep,pattern="YEAR")==1, arr.ind=TRUE)
## print PDO metadata
print(all_NPGO[seq(hdr_NPGO)],quote=FALSE)
##  [1]                                                                                                     
##  [2] <html>                                                                                              
##  [3] <body>                                                                                              
##  [4]                                                                                                     
##  [5] <pre>#  Last update 29-Mar-2017 by E. Di Lorenzo                                                    
##  [6] #  NPGO index monthly averages                                                                      
##  [7] #  from Jan-1950  to  Feb-2017                                                                      
##  [8] #                                                                                                   
##  [9] #  WARNING: Values after Dec-2004 are updated                                                       
## [10] #  using Satellite SSHa from AVISO Delayed Time product.                                            
## [11] #  http://opendap.aviso.oceanobs.com/thredds/dodsC/dataset-duacs-dt-global-allsat-msla-h            
## [12] #                                                                                                   
## [13] #  PRELIMINARY: Values after Jan-2016 are preliminary and updated                                   
## [14] #  using Satellite SSHa from AVISO Near Real Time product.                                          
## [15] #  http://opendap.aviso.oceanobs.com/thredds/dodsC/dataset-duacs-nrt-over30d-global-allsat-msla-h   
## [16] #                                                                                                   
## [17] #  The update is performed by taking the NPGO spatial pattern of Di Lorenzo et al. 2008             
## [18] #  computed over the period 1950-2004, and projecting the AVISO Satellite SSHa.                     
## [19] #  During the pre-processing of the AVISO data, we remove the seasonal cycle based on               
## [20] #  the 1993-2004 seasonal means.                                                                    
## [21] #                                                                                                   
## [22] #  AVISO PRODUCT UPDATE Summer 2014: AVISO has released a re-processed dataset for the sea level.   
## [23] #  Starting from the November 2014, the NPGO index is computed with this updated dataset. NPGO      
## [24] #  values from 2004 onward have been recomputed with very minor differences from previous releases. 
## [25] #                                                                                                   
## [26] #  Ref:                                                                                             
## [27] #  Di Lorenzo et al., 2008: North Pacific Gyre Oscillation                                          
## [28] #  links ocean climate and ecosystem change, GRL.                                                   
## [29] #                                                                                                   
## [30] #      YEAR            MONTH        NPGO index

Next, we will extract the actual NPGO indices for the years of interest and inspect the file contents. We also want the average annual NPGO index from January 1 through December 31 during the first year that the juvenile salmon are in the ocean (i.e., during their second year of life). Therefore, we need NPGO values from yr_first + 2 = 1962 through yr_last + 2 = 2013.

## NPGO data for years of interest
raw_NPGO <- read_table(url_NPGO, col_names = FALSE, skip=hdr_NPGO + (yr_first+2-1950)*12, n_max = TT*12)
colnames(raw_NPGO) <- c("year","month","NPGO")
## calculate annual means
dat_NPGO <- aggregate(raw_NPGO$NPGO, by = list(year = raw_NPGO$year), FUN = mean)
colnames(dat_NPGO) <- c("brood_yr", "NPGO")
## match calendar yr to brood year
dat_NPGO$brood_yr <- dat_NPGO$brood_yr - 2
dat_NPGO
## market index
market <- round(dat_NPGO$NPGO, 3)
plot.ts(market)

cor(market, t(assets), use="pairwise.complete.obs")
##    brood_yr         NPGO
## 1      1960 -0.093424126
## 2      1961 -0.234721502
## 3      1962  0.544102425
## 4      1963 -0.837212701
## 5      1964 -0.956356193
## 6      1965 -0.699724801
## 7      1966 -0.850601763
## 8      1967 -0.387628509
## 9      1968  0.061043154
## 10     1969  0.381984271
## 11     1970 -0.659303726
## 12     1971  0.126585564
## 13     1972 -0.009310217
## 14     1973  0.798500537
## 15     1974  1.863211350
## 16     1975  1.249270341
## 17     1976  0.504466636
## 18     1977 -0.737188575
## 19     1978 -0.749107276
## 20     1979 -0.221840562
## 21     1980 -0.686691381
## 22     1981 -0.092064694
## 23     1982  0.620750037
## 24     1983 -0.393901047
## 25     1984 -0.732847288
## 26     1985  0.429174729
## 27     1986  1.323636763
## 28     1987  0.612240617
## 29     1988 -0.201546121
## 30     1989 -0.498914711
## 31     1990 -1.094918870
## 32     1991 -1.901992588
## 33     1992 -1.211598783
## 34     1993 -1.039208546
## 35     1994 -0.972836612
## 36     1995 -0.657690616
## 37     1996  0.544553283
## 38     1997  1.515780517
## 39     1998  1.885861292
## 40     1999  2.080693375
## 41     2000  1.473593960
## 42     2001  0.958507189
## 43     2002  0.209460249
## 44     2003 -1.381917947
## 45     2004 -0.555041996
## 46     2005  0.593098207
## 47     2006  1.455274563
## 48     2007  0.736914838
##          GR_IR       GRCAT     GRLOS     GRMIN     GRUMA      GRWEN     MFBEA     MFBIG     MFCAM
## [1,] 0.2112462 -0.03420949 0.1016297 0.1658748 0.1063399 0.07647642 0.3068064 0.3254104 0.3094061
##          MFLOO     MFMAR     MFSUL     SFEFS     SFMAI     SFSEC     SREFS     SRLEM     SRLMA
## [1,] 0.3418419 0.1306568 0.2288056 0.2960428 0.2557321 0.4550384 0.3679295 0.1945687 0.3956831
##          SRUMA     SRVAL     SRYFS
## [1,] 0.4925672 0.3156856 0.3894858

2.3 Risk-free index

For these purposes, the risk free index is zero for all years. The idea is that any population with a standardized return of zero is replacing itself and therefore not declining.

free <- rep(0, TT)

3 CAPM in MARSS

3.1 v1: a static; b unconstrained

## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL

## number of regr params = level + slope
mm <- 2

## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
  
  ## print loop ID
  print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
  ## get popns
  assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
  ## number of assets
  nn <- dim(assets_mpg)[1]
  
  ## Process eqn
  ## B is Identity
  BB <- "identity"
  ## U is col vec of 0's
  UU <- "zero"
  ## Q will be block diagonal
  QQ <- matrix(list(0),mm*nn,mm*nn)
  ## upper L block for alpha
  aa <- 1:nn
  ## lower R block for beta
  bb <- aa + nn
  ## equal var-cov
  # QQ[bb,bb] <- "beta.cov"
  # diag(QQ[bb,bb]) <- "beta.var"
  ## diagonal & unequal
  diag(QQ[bb,bb]) <- paste0("beta",seq(nn))
  for(j in 1:(nn-1)) {
    for(k in (j+1):nn) {
      QQ[j+nn,k+nn] <- paste0("beta",j,k)
      QQ[k+nn,j+nn] <- paste0("beta",j,k)
    }
  }
  
  ## Observation eqn
  ## Z is array of regr variables (ie, 1's and market returns)
  ZZ <- array(NA, c(nn,mm*nn,TT))
  for(t in 1:TT) {
    ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
  }
  ## A is col vec of 0's
  AA <- "zero"
  ## assume same obs var & cov
  ##RR <- matrix(list(0),nn,nn)
  ##RR[aa,aa] <- "obs.cov"
  ##diag(RR[aa,aa]) <- "obs.var"
  RR <- "diagonal and equal"
  
  ## Initial states
  ##x0 <- rbind(matrix("pi.a",nn,1), matrix("pi.b",nn,1))
  ##V0 <- 100*diag(mm*nn)
  ##x0 <- rbind(matrix(0,nn,1), matrix(1,nn,1))
  ##V0 <- matrix(list(0),mm*nn,mm*nn)
  ##diag(V0[aa,aa]) <- "var.a"
  ##diag(V0[bb,bb]) <- "var.b"
  
  ## starting values for regr parameters
  ini_list <- list(x0=matrix(1, mm*nn, 1))
  
  ## list of model matrices & vectors
  mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
  ##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
  
  ## list of control values
  con_list <- list(safe=TRUE, allow.degen=FALSE,
                   conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
  
  ## fit multivariate DLM
  mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
  
  ## get list of Kalman filter output
  kf_out = MARSSkfss(mod_res[[i]])
  
  ## ts of forecasted betas
  betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
  
  ## ts of SE of betas; nxT matrix
  betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
  
  ## ts of alphas
  alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
  ## ts of SE{alphas}
  alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
  
  ## ts of betas
  betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
  ## ts of SE{betas}
  betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
  
}
## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 1500 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1500 iterations. 
## Log-likelihood: -390.0471 
## AIC: 848.0942   AICc: 857.8085   
##  
##          Estimate
## R.diag      0.706
## Q.beta1     0.719
## Q.beta12    0.942
## Q.beta13    0.756
## Q.beta14    0.604
## Q.beta15    0.872
## Q.beta16    0.762
## Q.beta2     1.234
## Q.beta23    0.991
## Q.beta24    0.791
## Q.beta25    1.142
## Q.beta26    0.999
## Q.beta3     0.796
## Q.beta34    0.635
## Q.beta35    0.917
## Q.beta36    0.802
## Q.beta4     0.507
## Q.beta45    0.732
## Q.beta46    0.640
## Q.beta5     1.057
## Q.beta56    0.924
## Q.beta6     0.808
## x0.X1      -0.287
## x0.X2      -0.326
## x0.X3      -0.158
## x0.X4      -0.174
## x0.X5      -0.344
## x0.X6      -0.293
## x0.X7      -0.743
## x0.X8      -1.343
## x0.X9      -0.880
## x0.X10     -0.614
## x0.X11     -1.034
## x0.X12     -0.943
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 1374 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1374 iterations. 
## Log-likelihood: -397.2663 
## AIC: 862.5327   AICc: 872.7036   
##  
##          Estimate
## R.diag      0.623
## Q.beta1     5.100
## Q.beta12    5.551
## Q.beta13    3.435
## Q.beta14    3.424
## Q.beta15    5.526
## Q.beta16    4.243
## Q.beta2     6.042
## Q.beta23    3.743
## Q.beta24    3.731
## Q.beta25    6.014
## Q.beta26    4.616
## Q.beta3     2.394
## Q.beta34    2.379
## Q.beta35    3.707
## Q.beta36    2.830
## Q.beta4     2.364
## Q.beta45    3.697
## Q.beta46    2.824
## Q.beta5     5.990
## Q.beta56    4.602
## Q.beta6     3.538
## x0.X1      -0.333
## x0.X2      -0.374
## x0.X3      -0.451
## x0.X4      -0.403
## x0.X5      -0.358
## x0.X6      -0.365
## x0.X7      -0.358
## x0.X8      -0.333
## x0.X9      -0.288
## x0.X10     -0.230
## x0.X11     -0.446
## x0.X12     -0.294
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 1013 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1013 iterations. 
## Log-likelihood: -167.7804 
## AIC: 361.5608   AICc: 364.5444   
##  
##          Estimate
## R.diag      0.370
## Q.beta1     1.405
## Q.beta12    1.076
## Q.beta13    1.112
## Q.beta2     0.825
## Q.beta23    0.852
## Q.beta3     0.881
## x0.X1      -0.277
## x0.X2      -0.237
## x0.X3      -0.143
## x0.X4       0.825
## x0.X5       0.636
## x0.X6       0.825
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 1276 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1276 iterations. 
## Log-likelihood: -371.8351 
## AIC: 811.6702   AICc: 821.4645   
##  
##          Estimate
## R.diag      0.381
## Q.beta1     9.543
## Q.beta12    7.379
## Q.beta13    6.356
## Q.beta14    7.864
## Q.beta15    8.931
## Q.beta16   12.867
## Q.beta2     5.796
## Q.beta23    4.951
## Q.beta24    6.078
## Q.beta25    6.838
## Q.beta26    9.838
## Q.beta3     4.249
## Q.beta34    5.238
## Q.beta35    5.923
## Q.beta36    8.528
## Q.beta4     6.482
## Q.beta45    7.364
## Q.beta46   10.611
## Q.beta5     8.413
## Q.beta56   12.131
## Q.beta6    17.495
## x0.X1      -0.645
## x0.X2      -0.696
## x0.X3      -0.566
## x0.X4      -0.535
## x0.X5      -0.576
## x0.X6      -0.779
## x0.X7      -5.837
## x0.X8      -4.579
## x0.X9      -3.728
## x0.X10     -4.512
## x0.X11     -5.483
## x0.X12     -7.824
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop

3.1.1 Plots of results

betas <- betas[,names_mat$ID]

# total num of ts
NN <- dim(betas)[2]

# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
             rainbow(3, start=0.2, end=0.45, v=0.7),
             rainbow(6, start=0.8, end=0.95, v=0.85),
             rainbow(6, start=0.55, end=0.7))

# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7

# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))

# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)

# draw plots
for(i in 1:NN) {
    plot(t_index, betas[,i], type="n",
         ylim=yext, xaxt="n", yaxt="n",
         xlab="", ylab="", main=names_mat[i,"long"])
    rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
    abline(h=0, lty="dashed")
    box()
    lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
    lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
    lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
    #text(xin, yin, names.pop[i], adj=c(0,1))   
    if(i<=nr | nc==1) {
#       axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
#       axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
        axis(side=2, las=1, tick=TRUE)
        }
    if(i%%nr==0) {
        axis(side=1, at=t_index[(t_index %% 10)==0])
        }
    }
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)

3.2 v1: a static; b equalvarcov

## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL

## number of regr params = level + slope
mm <- 2

## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
  
  ## print loop ID
  print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
  ## get popns
  assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
  ## number of assets
  nn <- dim(assets_mpg)[1]
  
  ## Process eqn
  ## B is Identity
  BB <- "identity"
  ## U is col vec of 0's
  UU <- "zero"
  ## Q will be block diagonal
  QQ <- matrix(list(0),mm*nn,mm*nn)
  ## upper L block for alpha
  aa <- 1:nn
  ## lower R block for beta
  bb <- aa + nn
  ## equal var-cov
  QQ[bb,bb] <- "beta.cov"
  diag(QQ[bb,bb]) <- "beta.var"

  ## Observation eqn
  ## Z is array of regr variables (ie, 1's and market returns)
  ZZ <- array(NA, c(nn,mm*nn,TT))
  for(t in 1:TT) {
    ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
  }
  ## A is col vec of 0's
  AA <- "zero"
  ## assume same obs var & cov
  ##RR <- matrix(list(0),nn,nn)
  ##RR[aa,aa] <- "obs.cov"
  ##diag(RR[aa,aa]) <- "obs.var"
  RR <- "diagonal and equal"
  
  ## Initial states
  ##x0 <- rbind(matrix("pi.a",nn,1), matrix("pi.b",nn,1))
  ##V0 <- 100*diag(mm*nn)
  ##x0 <- rbind(matrix(0,nn,1), matrix(1,nn,1))
  ##V0 <- matrix(list(0),mm*nn,mm*nn)
  ##diag(V0[aa,aa]) <- "var.a"
  ##diag(V0[bb,bb]) <- "var.b"
  
  ## starting values for regr parameters
  ini_list <- list(x0=matrix(1, mm*nn, 1))
  
  ## list of model matrices & vectors
  mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
  ##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
  
  ## list of control values
  con_list <- list(safe=TRUE, allow.degen=FALSE,
                   conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
  
  ## fit multivariate DLM
  mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
  
  ## get list of Kalman filter output
  kf_out = MARSSkfss(mod_res[[i]])
  
  ## ts of forecasted betas
  betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
  
  ## ts of SE of betas; nxT matrix
  betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
  
  ## ts of alphas
  alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
  ## ts of SE{alphas}
  alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
  
  ## ts of betas
  betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
  ## ts of SE{betas}
  betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
  
}
## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 1415 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1415 iterations. 
## Log-likelihood: -392.2293 
## AIC: 814.4587   AICc: 816.2769   
##  
##            Estimate
## R.diag        0.720
## Q.beta.var    0.833
## Q.beta.cov    0.833
## x0.X1        -0.296
## x0.X2        -0.269
## x0.X3        -0.157
## x0.X4        -0.211
## x0.X5        -0.315
## x0.X6        -0.293
## x0.X7        -0.879
## x0.X8        -1.175
## x0.X9        -0.969
## x0.X10       -0.903
## x0.X11       -0.967
## x0.X12       -1.022
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 1338 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1338 iterations. 
## Log-likelihood: -399.7483 
## AIC: 829.4966   AICc: 831.3938   
##  
##            Estimate
## R.diag       0.7063
## Q.beta.var   3.6593
## Q.beta.cov   3.6593
## x0.X1       -0.3051
## x0.X2       -0.3234
## x0.X3       -0.4413
## x0.X4       -0.4151
## x0.X5       -0.3175
## x0.X6       -0.3518
## x0.X7       -0.1535
## x0.X8       -0.0502
## x0.X9       -0.1131
## x0.X10      -0.0547
## x0.X11      -0.1856
## x0.X12      -0.2453
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 984 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 984 iterations. 
## Log-likelihood: -169.5842 
## AIC: 357.1684   AICc: 358.597   
##  
##            Estimate
## R.diag        0.397
## Q.beta.var    0.908
## Q.beta.cov    0.908
## x0.X1        -0.246
## x0.X2        -0.239
## x0.X3        -0.141
## x0.X4         0.712
## x0.X5         0.645
## x0.X6         0.818
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 1458 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1458 iterations. 
## Log-likelihood: -387.8977 
## AIC: 805.7954   AICc: 807.6275   
##  
##            Estimate
## R.diag        0.522
## Q.beta.var    6.799
## Q.beta.cov    6.799
## x0.X1        -0.505
## x0.X2        -0.611
## x0.X3        -0.600
## x0.X4        -0.507
## x0.X5        -0.549
## x0.X6        -0.635
## x0.X7        -3.115
## x0.X8        -3.437
## x0.X9        -3.266
## x0.X10       -3.025
## x0.X11       -3.189
## x0.X12       -2.944
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop

3.2.1 Plots of results

betas <- betas[,names_mat$ID]

# total num of ts
NN <- dim(betas)[2]

# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
             rainbow(3, start=0.2, end=0.45, v=0.7),
             rainbow(6, start=0.8, end=0.95, v=0.85),
             rainbow(6, start=0.55, end=0.7))

# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7

# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))

# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)

# draw plots
for(i in 1:NN) {
    plot(t_index, betas[,i], type="n",
         ylim=yext, xaxt="n", yaxt="n",
         xlab="", ylab="", main=names_mat[i,"long"])
    rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
    abline(h=0, lty="dashed")
    box()
    lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
    lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
    lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
    #text(xin, yin, names.pop[i], adj=c(0,1))   
    if(i<=nr | nc==1) {
#       axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
#       axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
        axis(side=2, las=1, tick=TRUE)
        }
    if(i%%nr==0) {
        axis(side=1, at=t_index[(t_index %% 10)==0])
        }
    }
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)

3.3 v1: a static; b diag & unequal

## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL

## number of regr params = level + slope
mm <- 2

## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
  
  ## print loop ID
  print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
  ## get popns
  assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
  ## number of assets
  nn <- dim(assets_mpg)[1]
  
  ## Process eqn
  ## B is Identity
  BB <- "identity"
  ## U is col vec of 0's
  UU <- "zero"
  ## Q will be block diagonal
  QQ <- matrix(list(0),mm*nn,mm*nn)
  ## upper L block for alpha
  aa <- 1:nn
  ## lower R block for beta
  bb <- aa + nn
  ## diag & unequal
  diag(QQ[bb,bb]) <- paste0("beta",seq(nn))

  ## Observation eqn
  ## Z is array of regr variables (ie, 1's and market returns)
  ZZ <- array(NA, c(nn,mm*nn,TT))
  for(t in 1:TT) {
    ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
  }
  ## A is col vec of 0's
  AA <- "zero"
  ## assume same obs var & cov
  ##RR <- matrix(list(0),nn,nn)
  ##RR[aa,aa] <- "obs.cov"
  ##diag(RR[aa,aa]) <- "obs.var"
  RR <- "diagonal and equal"
  
  ## Initial states
  ##x0 <- rbind(matrix("pi.a",nn,1), matrix("pi.b",nn,1))
  ##V0 <- 100*diag(mm*nn)
  ##x0 <- rbind(matrix(0,nn,1), matrix(1,nn,1))
  ##V0 <- matrix(list(0),mm*nn,mm*nn)
  ##diag(V0[aa,aa]) <- "var.a"
  ##diag(V0[bb,bb]) <- "var.b"
  
  ## starting values for regr parameters
  ini_list <- list(x0=matrix(1, mm*nn, 1))
  
  ## list of model matrices & vectors
  mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
  ##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
  
  ## list of control values
  con_list <- list(safe=TRUE, allow.degen=FALSE,
                   conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
  
  ## fit multivariate DLM
  mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
  
  ## get list of Kalman filter output
  kf_out = MARSSkfss(mod_res[[i]])
  
  ## ts of forecasted betas
  betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
  
  ## ts of SE of betas; nxT matrix
  betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
  
  ## ts of alphas
  alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
  ## ts of SE{alphas}
  alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
  
  ## ts of betas
  betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
  ## ts of SE{betas}
  betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
  
}
## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 1534 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1534 iterations. 
## Log-likelihood: -425.1244 
## AIC: 888.2489   AICc: 891.172   
##  
##          Estimate
## R.diag   1.11e+00
## Q.beta1  1.56e-02
## Q.beta2  8.14e-02
## Q.beta3  1.78e-02
## Q.beta4  4.08e-05
## Q.beta5  2.87e-02
## Q.beta6  1.08e-04
## x0.X1   -1.94e-01
## x0.X2   -2.68e-01
## x0.X3   -3.38e-02
## x0.X4   -2.97e-02
## x0.X5   -2.33e-01
## x0.X6   -1.18e-01
## x0.X7   -2.51e-01
## x0.X8   -1.19e+00
## x0.X9   -3.06e-02
## x0.X10   1.80e-01
## x0.X11  -7.18e-01
## x0.X12   7.81e-02
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 3014 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 3014 iterations. 
## Log-likelihood: -452.207 
## AIC: 942.414   AICc: 945.4662   
##  
##          Estimate
## R.diag   1.65e+00
## Q.beta1  4.15e-05
## Q.beta2  3.08e-03
## Q.beta3  2.30e-02
## Q.beta4  1.46e-04
## Q.beta5  5.88e-05
## Q.beta6  4.99e-05
## x0.X1   -1.13e-01
## x0.X2   -1.49e-01
## x0.X3   -3.36e-01
## x0.X4   -1.21e-01
## x0.X5   -1.53e-01
## x0.X6   -3.63e-01
## x0.X7    4.13e-01
## x0.X8    3.58e-01
## x0.X9   -7.52e-02
## x0.X10   5.15e-01
## x0.X11   1.84e-01
## x0.X12   3.57e-01
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 2082 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2082 iterations. 
## Log-likelihood: -181.13 
## AIC: 382.2601   AICc: 384.0201   
##  
##          Estimate
## R.diag   0.817099
## Q.beta1  0.008137
## Q.beta2  0.000028
## Q.beta3  0.000066
## x0.X1   -0.166970
## x0.X2   -0.120705
## x0.X3    0.005672
## x0.X4   -0.048574
## x0.X5    0.229948
## x0.X6    0.485909
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 1740 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1740 iterations. 
## Log-likelihood: -463.9245 
## AIC: 965.849   AICc: 968.7947   
##  
##          Estimate
## R.diag   1.65e+00
## Q.beta1  1.95e-04
## Q.beta2  4.43e-05
## Q.beta3  4.26e-05
## Q.beta4  4.39e-05
## Q.beta5  7.75e-05
## Q.beta6  1.89e-04
## x0.X1   -1.43e-01
## x0.X2   -2.48e-01
## x0.X3   -1.86e-01
## x0.X4   -1.44e-01
## x0.X5   -2.04e-01
## x0.X6   -3.86e-01
## x0.X7    5.58e-01
## x0.X8    2.43e-01
## x0.X9    4.70e-01
## x0.X10   6.56e-01
## x0.X11   4.73e-01
## x0.X12   7.64e-01
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop

3.3.1 Plots of results

betas <- betas[,names_mat$ID]

# total num of ts
NN <- dim(betas)[2]

# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
             rainbow(3, start=0.2, end=0.45, v=0.7),
             rainbow(6, start=0.8, end=0.95, v=0.85),
             rainbow(6, start=0.55, end=0.7))

# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7

# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))

# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)

# draw plots
for(i in 1:NN) {
    plot(t_index, betas[,i], type="n",
         ylim=yext, xaxt="n", yaxt="n",
         xlab="", ylab="", main=names_mat[i,"long"])
    rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
    abline(h=0, lty="dashed")
    box()
    lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
    lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
    lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
    #text(xin, yin, names.pop[i], adj=c(0,1))   
    if(i<=nr | nc==1) {
#       axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
#       axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
        axis(side=2, las=1, tick=TRUE)
        }
    if(i%%nr==0) {
        axis(side=1, at=t_index[(t_index %% 10)==0])
        }
    }
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)

3.4 v2: a unconstrained; b unconstrained

## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL

## number of regr params = level + slope
mm <- 2

## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
  
  ## print loop ID
  print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
  ## get popns
  assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
  ## number of assets
  nn <- dim(assets_mpg)[1]
  
  ## Process eqn
  ## B is Identity
  BB <- "identity"
  ## U is col vec of 0's
  UU <- "zero"
  ## Q will be block diagonal
  QQ <- matrix(list(0),mm*nn,mm*nn)
  ## upper L block for alpha
  aa <- 1:nn
  ## equal var-cov
  # QQ[aa,aa] <- "alpha.cov"
  # diag(QQ[aa,aa]) <- "alpha.var"
  ## diagonal & unequal
  diag(QQ[aa,aa]) <- paste0("alpha",seq(nn))
  for(j in 1:(nn-1)) {
    for(k in (j+1):nn) {
      QQ[j,k] <- paste0("alpha",j,k)
      QQ[k,j] <- paste0("alpha",j,k)
    }
  }
  ## lower R block for beta
  bb <- aa + nn
  ## equal var-cov
  # QQ[bb,bb] <- "beta.cov"
  # diag(QQ[bb,bb]) <- "beta.var"
  ## diagonal & unequal
  diag(QQ[bb,bb]) <- paste0("beta",seq(nn))
  for(j in 1:(nn-1)) {
    for(k in (j+1):nn) {
      QQ[j+nn,k+nn] <- paste0("beta",j,k)
      QQ[k+nn,j+nn] <- paste0("beta",j,k)
    }
  }
  
  ## Observation eqn
  ## Z is array of regr variables (ie, 1's and market returns)
  ZZ <- array(NA, c(nn,mm*nn,TT))
  for(t in 1:TT) {
    ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
  }
  ## A is col vec of 0's
  AA <- "zero"
  ## assume same obs var & cov
  ##RR <- matrix(list(0),nn,nn)
  ##RR[aa,aa] <- "obs.cov"
  ##diag(RR[aa,aa]) <- "obs.var"
  RR <- "diagonal and equal"
  
  ## starting values for regr parameters
  ini_list <- list(x0=matrix(1, mm*nn, 1))
  
  ## list of model matrices & vectors
  mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
  ##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
  
  ## list of control values
  con_list <- list(safe=TRUE, allow.degen=FALSE,
                   conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
  
  ## fit multivariate DLM
  mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
  
  ## get list of Kalman filter output
  kf_out = MARSSkfss(mod_res[[i]])
  
  ## ts of forecasted betas
  betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
  
  ## ts of SE of betas; nxT matrix
  betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
  
  ## ts of alphas
  alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
  ## ts of SE{alphas}
  alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
  
  ## ts of betas
  betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
  ## ts of SE{betas}
  betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
  
}
## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 2703 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2703 iterations. 
## Log-likelihood: -363.2693 
## AIC: 836.5386   AICc: 864.0386   
##  
##           Estimate
## R.diag      0.3003
## Q.alpha1    0.4947
## Q.alpha12   0.5072
## Q.alpha13   0.6028
## Q.alpha14   0.4998
## Q.alpha15   0.2693
## Q.alpha16   0.4113
## Q.alpha2    0.7917
## Q.alpha23   0.5354
## Q.alpha24   0.3092
## Q.alpha25   0.3237
## Q.alpha26   0.3896
## Q.alpha3    0.7968
## Q.alpha34   0.6822
## Q.alpha35   0.4049
## Q.alpha36   0.5219
## Q.alpha4    0.6614
## Q.alpha45   0.2649
## Q.alpha46   0.4461
## Q.alpha5    0.3794
## Q.alpha56   0.2462
## Q.alpha6    0.3588
## Q.beta1     0.0132
## Q.beta12    0.0388
## Q.beta13    0.0140
## Q.beta14    0.0324
## Q.beta15    0.0979
## Q.beta16    0.0535
## Q.beta2     0.1244
## Q.beta23    0.0465
## Q.beta24    0.1163
## Q.beta25    0.3303
## Q.beta26    0.1848
## Q.beta3     0.0177
## Q.beta34    0.0454
## Q.beta35    0.1260
## Q.beta36    0.0711
## Q.beta4     0.1226
## Q.beta45    0.3271
## Q.beta46    0.1876
## Q.beta5     0.9014
## Q.beta56    0.5103
## Q.beta6     0.2903
## x0.X1      -0.4679
## x0.X2      -1.7161
## x0.X3       0.3303
## x0.X4       0.5274
## x0.X5       0.0813
## x0.X6      -0.5697
## x0.X7      -0.1207
## x0.X8      -0.9706
## x0.X9      -0.0274
## x0.X10     -0.0442
## x0.X11     -1.5186
## x0.X12     -0.6445
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 2135 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2135 iterations. 
## Log-likelihood: -352.7901 
## AIC: 815.5802   AICc: 844.5004   
##  
##            Estimate
## R.diag     0.281363
## Q.alpha1   1.429407
## Q.alpha12  1.613941
## Q.alpha13  0.673732
## Q.alpha14  0.692343
## Q.alpha15  1.489888
## Q.alpha16  1.608054
## Q.alpha2   1.842574
## Q.alpha23  0.797055
## Q.alpha24  0.823314
## Q.alpha25  1.681005
## Q.alpha26  1.776505
## Q.alpha3   0.924359
## Q.alpha34  1.097678
## Q.alpha35  0.720479
## Q.alpha36  0.109656
## Q.alpha4   1.317354
## Q.alpha45  0.745323
## Q.alpha46 -0.045377
## Q.alpha5   1.553820
## Q.alpha56  1.656737
## Q.alpha6   2.501869
## Q.beta1    0.002804
## Q.beta12   0.001372
## Q.beta13  -0.002389
## Q.beta14   0.002082
## Q.beta15   0.006010
## Q.beta16  -0.006659
## Q.beta2    0.000708
## Q.beta23  -0.001106
## Q.beta24   0.001024
## Q.beta25   0.002905
## Q.beta26  -0.003139
## Q.beta3    0.002238
## Q.beta34  -0.001723
## Q.beta35  -0.005229
## Q.beta36   0.005939
## Q.beta4    0.001635
## Q.beta45   0.004485
## Q.beta46  -0.005042
## Q.beta5    0.012994
## Q.beta56  -0.014538
## Q.beta6    0.016617
## x0.X1      0.231601
## x0.X2     -0.288871
## x0.X3      0.008899
## x0.X4      0.223452
## x0.X5      0.335920
## x0.X6      0.100134
## x0.X7      0.509994
## x0.X8      0.522791
## x0.X9      0.094360
## x0.X10     0.126260
## x0.X11     0.478596
## x0.X12     0.416539
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 2136 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2136 iterations. 
## Log-likelihood: -135.2344 
## AIC: 308.4688   AICc: 315.0206   
##  
##            Estimate
## R.diag     1.37e-01
## Q.alpha1   7.71e-01
## Q.alpha12  6.34e-01
## Q.alpha13  6.02e-01
## Q.alpha2   5.23e-01
## Q.alpha23  4.89e-01
## Q.alpha3   5.48e-01
## Q.beta1    1.80e-04
## Q.beta12   1.01e-04
## Q.beta13   1.06e-04
## Q.beta2    6.48e-05
## Q.beta23   6.18e-05
## Q.beta3    7.57e-05
## x0.X1     -9.75e-01
## x0.X2     -8.61e-01
## x0.X3     -6.27e-01
## x0.X4      3.00e-01
## x0.X5      2.41e-01
## x0.X6      4.25e-01
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 1577 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1577 iterations. 
## Log-likelihood: -338.4211 
## AIC: 786.8422   AICc: 814.59   
##  
##           Estimate
## R.diag     0.28630
## Q.alpha1   1.34857
## Q.alpha12  0.96660
## Q.alpha13  0.95088
## Q.alpha14  1.06089
## Q.alpha15  1.28876
## Q.alpha16  1.68622
## Q.alpha2   0.69572
## Q.alpha23  0.67684
## Q.alpha24  0.77004
## Q.alpha25  0.93189
## Q.alpha26  1.22611
## Q.alpha3   0.67868
## Q.alpha34  0.73401
## Q.alpha35  0.89737
## Q.alpha36  1.16775
## Q.alpha4   0.87348
## Q.alpha45  1.04886
## Q.alpha46  1.41327
## Q.alpha5   1.26369
## Q.alpha56  1.69375
## Q.alpha6   2.33465
## Q.beta1    0.00653
## Q.beta12   0.03641
## Q.beta13   0.01717
## Q.beta14   0.00838
## Q.beta15  -0.01926
## Q.beta16  -0.01446
## Q.beta2    0.20737
## Q.beta23   0.09574
## Q.beta24   0.04554
## Q.beta25  -0.11392
## Q.beta26  -0.08513
## Q.beta3    0.04585
## Q.beta34   0.02275
## Q.beta35  -0.04938
## Q.beta36  -0.03737
## Q.beta4    0.01184
## Q.beta45  -0.02161
## Q.beta46  -0.01663
## Q.beta5    0.06896
## Q.beta56   0.05064
## Q.beta6    0.03738
## x0.X1     -0.08959
## x0.X2     -0.17092
## x0.X3     -0.27497
## x0.X4      0.36118
## x0.X5      0.38791
## x0.X6      1.08623
## x0.X7      0.49013
## x0.X8      0.10929
## x0.X9      0.49369
## x0.X10     0.71621
## x0.X11     0.66461
## x0.X12     0.88863
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop

3.4.1 Plots of results

betas <- betas[,names_mat$ID]

# total num of ts
NN <- dim(betas)[2]

# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
             rainbow(3, start=0.2, end=0.45, v=0.7),
             rainbow(6, start=0.8, end=0.95, v=0.85),
             rainbow(6, start=0.55, end=0.7))

# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7

# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))

# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)

# draw plots
for(i in 1:NN) {
    plot(t_index, betas[,i], type="n",
         ylim=yext, xaxt="n", yaxt="n",
         xlab="", ylab="", main=names_mat[i,"long"])
    rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
    abline(h=0, lty="dashed")
    box()
    lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
    lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
    lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
    #text(xin, yin, names.pop[i], adj=c(0,1))   
    if(i<=nr | nc==1) {
#       axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
#       axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
        axis(side=2, las=1, tick=TRUE)
        }
    if(i%%nr==0) {
        axis(side=1, at=t_index[(t_index %% 10)==0])
        }
    }
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)

3.5 v2: a unconstrained; b equalvarcov

## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL

## number of regr params = level + slope
mm <- 2

## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
  
  ## print loop ID
  print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
  ## get popns
  assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
  ## number of assets
  nn <- dim(assets_mpg)[1]
  
  ## Process eqn
  ## B is Identity
  BB <- "identity"
  ## U is col vec of 0's
  UU <- "zero"
  ## Q will be block diagonal
  QQ <- matrix(list(0),mm*nn,mm*nn)
  ## upper L block for alpha
  aa <- 1:nn
  ## equal var-cov
  # QQ[aa,aa] <- "alpha.cov"
  # diag(QQ[aa,aa]) <- "alpha.var"
  ## diagonal & unequal
  diag(QQ[aa,aa]) <- paste0("alpha",seq(nn))
  for(j in 1:(nn-1)) {
    for(k in (j+1):nn) {
      QQ[j,k] <- paste0("alpha",j,k)
      QQ[k,j] <- paste0("alpha",j,k)
    }
  }
  ## lower R block for beta
  bb <- aa + nn
  ## equal var-cov
  QQ[bb,bb] <- "beta.cov"
  diag(QQ[bb,bb]) <- "beta.var"
  ## diagonal & unequal
  # diag(QQ[bb,bb]) <- paste0("beta",seq(nn))
  # for(j in 1:(nn-1)) {
  #   for(k in (j+1):nn) {
  #     QQ[j+nn,k+nn] <- paste0("beta",j,k)
  #     QQ[k+nn,j+nn] <- paste0("beta",j,k)
  #   }
  # }
  
  ## Observation eqn
  ## Z is array of regr variables (ie, 1's and market returns)
  ZZ <- array(NA, c(nn,mm*nn,TT))
  for(t in 1:TT) {
    ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
  }
  ## A is col vec of 0's
  AA <- "zero"
  ## assume same obs var & cov
  ##RR <- matrix(list(0),nn,nn)
  ##RR[aa,aa] <- "obs.cov"
  ##diag(RR[aa,aa]) <- "obs.var"
  RR <- "diagonal and equal"
  
  ## starting values for regr parameters
  ini_list <- list(x0=matrix(1, mm*nn, 1))
  
  ## list of model matrices & vectors
  mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
  ##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
  
  ## list of control values
  con_list <- list(safe=TRUE, allow.degen=FALSE,
                   conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
  
  ## fit multivariate DLM
  mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
  
  ## get list of Kalman filter output
  kf_out = MARSSkfss(mod_res[[i]])
  
  ## ts of forecasted betas
  betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
  
  ## ts of SE of betas; nxT matrix
  betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
  
  ## ts of alphas
  alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
  ## ts of SE{alphas}
  alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
  
  ## ts of betas
  betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
  ## ts of SE{betas}
  betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
  
}
## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 1956 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1956 iterations. 
## Log-likelihood: -367.3308 
## AIC: 806.6615   AICc: 817.6245   
##  
##             Estimate
## R.diag      4.70e-01
## Q.alpha1    4.69e-01
## Q.alpha12   5.73e-01
## Q.alpha13   5.86e-01
## Q.alpha14   4.67e-01
## Q.alpha15   5.27e-01
## Q.alpha16   4.84e-01
## Q.alpha2    8.52e-01
## Q.alpha23   6.45e-01
## Q.alpha24   4.34e-01
## Q.alpha25   6.06e-01
## Q.alpha26   5.36e-01
## Q.alpha3    7.73e-01
## Q.alpha34   6.50e-01
## Q.alpha35   6.87e-01
## Q.alpha36   6.25e-01
## Q.alpha4    5.89e-01
## Q.alpha45   5.63e-01
## Q.alpha46   5.29e-01
## Q.alpha5    6.15e-01
## Q.alpha56   5.51e-01
## Q.alpha6    5.20e-01
## Q.beta.var  6.89e-05
## Q.beta.cov  5.03e-05
## x0.X1      -3.95e-01
## x0.X2      -1.24e+00
## x0.X3       2.67e-01
## x0.X4       5.05e-01
## x0.X5       2.21e-02
## x0.X6      -2.13e-01
## x0.X7       2.70e-01
## x0.X8      -1.02e-01
## x0.X9       2.08e-01
## x0.X10      2.88e-01
## x0.X11      1.81e-01
## x0.X12      1.64e-01
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 2051 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2051 iterations. 
## Log-likelihood: -353.0967 
## AIC: 778.1934   AICc: 789.6761   
##  
##             Estimate
## R.diag      3.06e-01
## Q.alpha1    1.44e+00
## Q.alpha12   1.59e+00
## Q.alpha13   6.86e-01
## Q.alpha14   7.06e-01
## Q.alpha15   1.50e+00
## Q.alpha16   1.63e+00
## Q.alpha2    1.78e+00
## Q.alpha23   8.02e-01
## Q.alpha24   8.34e-01
## Q.alpha25   1.66e+00
## Q.alpha26   1.76e+00
## Q.alpha3    9.56e-01
## Q.alpha34   1.11e+00
## Q.alpha35   7.30e-01
## Q.alpha36   1.43e-01
## Q.alpha4    1.29e+00
## Q.alpha45   7.54e-01
## Q.alpha46   2.23e-02
## Q.alpha5    1.57e+00
## Q.alpha56   1.69e+00
## Q.alpha6    2.49e+00
## Q.beta.var  6.00e-05
## Q.beta.cov  4.71e-05
## x0.X1       1.94e-01
## x0.X2      -2.56e-01
## x0.X3       1.09e-01
## x0.X4       1.72e-01
## x0.X5       2.56e-01
## x0.X6       1.56e-01
## x0.X7       4.48e-01
## x0.X8       4.99e-01
## x0.X9       1.64e-01
## x0.X10      1.01e-01
## x0.X11      3.59e-01
## x0.X12      5.38e-01
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 1564 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1564 iterations. 
## Log-likelihood: -135.2571 
## AIC: 300.5141   AICc: 304.5141   
##  
##             Estimate
## R.diag      1.37e-01
## Q.alpha1    7.77e-01
## Q.alpha12   6.35e-01
## Q.alpha13   6.04e-01
## Q.alpha2    5.19e-01
## Q.alpha23   4.87e-01
## Q.alpha3    5.48e-01
## Q.beta.var  4.78e-05
## Q.beta.cov  3.70e-05
## x0.X1      -9.73e-01
## x0.X2      -8.56e-01
## x0.X3      -6.27e-01
## x0.X4       3.09e-01
## x0.X5       2.44e-01
## x0.X6       4.31e-01
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 2028 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2028 iterations. 
## Log-likelihood: -340.9977 
## AIC: 753.9953   AICc: 765.0492   
##  
##             Estimate
## R.diag      3.56e-01
## Q.alpha1    1.32e+00
## Q.alpha12   1.00e+00
## Q.alpha13   9.35e-01
## Q.alpha14   1.07e+00
## Q.alpha15   1.23e+00
## Q.alpha16   1.67e+00
## Q.alpha2    7.79e-01
## Q.alpha23   7.13e-01
## Q.alpha24   7.91e-01
## Q.alpha25   8.92e-01
## Q.alpha26   1.20e+00
## Q.alpha3    6.63e-01
## Q.alpha34   7.52e-01
## Q.alpha35   8.63e-01
## Q.alpha36   1.17e+00
## Q.alpha4    8.80e-01
## Q.alpha45   1.03e+00
## Q.alpha46   1.42e+00
## Q.alpha5    1.23e+00
## Q.alpha56   1.71e+00
## Q.alpha6    2.39e+00
## Q.beta.var  8.35e-05
## Q.beta.cov  6.97e-05
## x0.X1       6.78e-03
## x0.X2      -4.22e-01
## x0.X3      -1.70e-01
## x0.X4       2.46e-01
## x0.X5       5.47e-01
## x0.X6       9.87e-01
## x0.X7       5.06e-01
## x0.X8       1.83e-01
## x0.X9       3.89e-01
## x0.X10      6.20e-01
## x0.X11      4.60e-01
## x0.X12      6.94e-01
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop

3.5.1 Plots of results

betas <- betas[,names_mat$ID]

# total num of ts
NN <- dim(betas)[2]

# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
             rainbow(3, start=0.2, end=0.45, v=0.7),
             rainbow(6, start=0.8, end=0.95, v=0.85),
             rainbow(6, start=0.55, end=0.7))

# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7

# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))

# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)

# draw plots
for(i in 1:NN) {
    plot(t_index, betas[,i], type="n",
         ylim=yext, xaxt="n", yaxt="n",
         xlab="", ylab="", main=names_mat[i,"long"])
    rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
    abline(h=0, lty="dashed")
    box()
    lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
    lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
    lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
    #text(xin, yin, names.pop[i], adj=c(0,1))   
    if(i<=nr | nc==1) {
#       axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
#       axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
        axis(side=2, las=1, tick=TRUE)
        }
    if(i%%nr==0) {
        axis(side=1, at=t_index[(t_index %% 10)==0])
        }
    }
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)

3.6 v3: a equalvarcov; b unconstrained

## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL

## number of regr params = level + slope
mm <- 2

## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
  
  ## print loop ID
  print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
  ## get popns
  assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
  ## number of assets
  nn <- dim(assets_mpg)[1]
  
  ## Process eqn
  ## B is Identity
  BB <- "identity"
  ## U is col vec of 0's
  UU <- "zero"
  ## Q will be block diagonal
  QQ <- matrix(list(0),mm*nn,mm*nn)
  ## upper L block for alpha
  aa <- 1:nn
  ## equal var-cov
  QQ[aa,aa] <- "alpha.cov"
  diag(QQ[aa,aa]) <- "alpha.var"
  ## lower R block for beta
  bb <- aa + nn
  ## diagonal & unequal
  diag(QQ[bb,bb]) <- paste0("beta",seq(nn))
  for(j in 1:(nn-1)) {
    for(k in (j+1):nn) {
      QQ[j+nn,k+nn] <- paste0("beta",j,k)
      QQ[k+nn,j+nn] <- paste0("beta",j,k)
    }
  }
  
  ## Observation eqn
  ## Z is array of regr variables (ie, 1's and market returns)
  ZZ <- array(NA, c(nn,mm*nn,TT))
  for(t in 1:TT) {
    ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
  }
  ## A is col vec of 0's
  AA <- "zero"
  ## assume same obs var & cov
  ##RR <- matrix(list(0),nn,nn)
  ##RR[aa,aa] <- "obs.cov"
  ##diag(RR[aa,aa]) <- "obs.var"
  RR <- "diagonal and equal"
  
  ## starting values for regr parameters
  ini_list <- list(x0=matrix(1, mm*nn, 1))
  
  ## list of model matrices & vectors
  mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
  ##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
  
  ## list of control values
  con_list <- list(safe=TRUE, allow.degen=FALSE,
                   conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
  
  ## fit multivariate DLM
  mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
  
  ## get list of Kalman filter output
  kf_out = MARSSkfss(mod_res[[i]])
  
  ## ts of forecasted betas
  betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
  
  ## ts of SE of betas; nxT matrix
  betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
  
  ## ts of alphas
  alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
  ## ts of SE{alphas}
  alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
  
  ## ts of betas
  betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
  ## ts of SE{betas}
  betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
  
}
## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 1962 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1962 iterations. 
## Log-likelihood: -371.5489 
## AIC: 815.0978   AICc: 826.0608   
##  
##              Estimate
## R.diag       0.512873
## Q.alpha.var  0.472812
## Q.alpha.cov  0.472788
## Q.beta1      0.012384
## Q.beta12     0.033173
## Q.beta13    -0.000533
## Q.beta14    -0.014850
## Q.beta15     0.005518
## Q.beta16    -0.001409
## Q.beta2      0.089998
## Q.beta23    -0.003417
## Q.beta24    -0.048698
## Q.beta25     0.001972
## Q.beta26    -0.012209
## Q.beta3      0.004044
## Q.beta34     0.018773
## Q.beta35     0.026132
## Q.beta36     0.017356
## Q.beta4      0.100114
## Q.beta45     0.113228
## Q.beta46     0.080249
## Q.beta5      0.177245
## Q.beta56     0.113874
## Q.beta6      0.075208
## x0.X1       -0.271026
## x0.X2       -0.350586
## x0.X3       -0.071601
## x0.X4       -0.085512
## x0.X5       -0.338093
## x0.X6       -0.222222
## x0.X7       -0.158583
## x0.X8       -1.184148
## x0.X9        0.055420
## x0.X10       0.015174
## x0.X11      -1.186145
## x0.X12      -0.568935
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 4799 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 4799 iterations. 
## Log-likelihood: -376.5283 
## AIC: 825.0567   AICc: 836.5394   
##  
##              Estimate
## R.diag       0.574195
## Q.alpha.var  1.090505
## Q.alpha.cov  1.090496
## Q.beta1      0.001277
## Q.beta12     0.002277
## Q.beta13     0.010294
## Q.beta14     0.009033
## Q.beta15    -0.000374
## Q.beta16    -0.005056
## Q.beta2      0.004094
## Q.beta23     0.018539
## Q.beta24     0.016265
## Q.beta25    -0.000690
## Q.beta26    -0.009128
## Q.beta3      0.084323
## Q.beta34     0.073969
## Q.beta35    -0.003216
## Q.beta36    -0.041622
## Q.beta4      0.064896
## Q.beta45    -0.002819
## Q.beta46    -0.036510
## Q.beta5      0.000152
## Q.beta56     0.001618
## Q.beta6      0.020596
## x0.X1        0.159421
## x0.X2        0.136142
## x0.X3       -0.000687
## x0.X4        0.044807
## x0.X5        0.158152
## x0.X6        0.105142
## x0.X7        0.433713
## x0.X8        0.531834
## x0.X9        0.377478
## x0.X10       0.415484
## x0.X11       0.342442
## x0.X12       0.325962
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 2022 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2022 iterations. 
## Log-likelihood: -139.7712 
## AIC: 309.5423   AICc: 313.5423   
##  
##             Estimate
## R.diag       0.20940
## Q.alpha.var  0.50907
## Q.alpha.cov  0.50906
## Q.beta1      0.00870
## Q.beta12     0.00481
## Q.beta13     0.00492
## Q.beta2      0.00267
## Q.beta23     0.00272
## Q.beta3      0.00279
## x0.X1       -0.90546
## x0.X2       -0.91729
## x0.X3       -0.78593
## x0.X4       -0.06701
## x0.X5        0.06959
## x0.X6        0.22662
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 1864 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1864 iterations. 
## Log-likelihood: -352.8772 
## AIC: 777.7545   AICc: 788.8084   
##  
##             Estimate
## R.diag       0.38467
## Q.alpha.var  0.97030
## Q.alpha.cov  0.97029
## Q.beta1      0.00124
## Q.beta12    -0.00892
## Q.beta13    -0.00240
## Q.beta14     0.00136
## Q.beta15     0.01100
## Q.beta16     0.02011
## Q.beta2      0.08153
## Q.beta23     0.01662
## Q.beta24    -0.02014
## Q.beta25    -0.11016
## Q.beta26    -0.20117
## Q.beta3      0.00545
## Q.beta34    -0.00109
## Q.beta35    -0.01852
## Q.beta36    -0.03411
## Q.beta4      0.00940
## Q.beta45     0.03298
## Q.beta46     0.05982
## Q.beta5      0.15644
## Q.beta56     0.28511
## Q.beta6      0.51972
## x0.X1        0.18842
## x0.X2        0.05143
## x0.X3        0.12064
## x0.X4        0.22757
## x0.X5        0.24395
## x0.X6        0.17975
## x0.X7        0.42536
## x0.X8        0.46287
## x0.X9        0.54532
## x0.X10       0.69797
## x0.X11       0.24589
## x0.X12       0.36009
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop

3.6.1 Plots of results

betas <- betas[,names_mat$ID]

# total num of ts
NN <- dim(betas)[2]

# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
             rainbow(3, start=0.2, end=0.45, v=0.7),
             rainbow(6, start=0.8, end=0.95, v=0.85),
             rainbow(6, start=0.55, end=0.7))

# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7

# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))

# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)

# draw plots
for(i in 1:NN) {
    plot(t_index, betas[,i], type="n",
         ylim=yext, xaxt="n", yaxt="n",
         xlab="", ylab="", main=names_mat[i,"long"])
    rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
    abline(h=0, lty="dashed")
    box()
    lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
    lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
    lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
    #text(xin, yin, names.pop[i], adj=c(0,1))   
    if(i<=nr | nc==1) {
#       axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
#       axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
        axis(side=2, las=1, tick=TRUE)
        }
    if(i%%nr==0) {
        axis(side=1, at=t_index[(t_index %% 10)==0])
        }
    }
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)

3.7 v3: a equalvarcov; b equalvarcov

## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL

## number of regr params = level + slope
mm <- 2

## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
  
  ## print loop ID
  print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
  ## get popns
  assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
  ## number of assets
  nn <- dim(assets_mpg)[1]
  
  ## Process eqn
  ## B is Identity
  BB <- "identity"
  ## U is col vec of 0's
  UU <- "zero"
  ## Q will be block diagonal
  QQ <- matrix(list(0),mm*nn,mm*nn)
  ## upper L block for alpha
  aa <- 1:nn
  ## equal var-cov
  QQ[aa,aa] <- "alpha.cov"
  diag(QQ[aa,aa]) <- "alpha.var"
  ## lower R block for beta
  bb <- aa + nn
  ## equal var-cov
  QQ[bb,bb] <- "beta.cov"
  diag(QQ[bb,bb]) <- "beta.var"

  ## Observation eqn
  ## Z is array of regr variables (ie, 1's and market returns)
  ZZ <- array(NA, c(nn,mm*nn,TT))
  for(t in 1:TT) {
    ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
  }
  ## A is col vec of 0's
  AA <- "zero"
  ## assume same obs var & cov
  ##RR <- matrix(list(0),nn,nn)
  ##RR[aa,aa] <- "obs.cov"
  ##diag(RR[aa,aa]) <- "obs.var"
  RR <- "diagonal and equal"
  
  ## starting values for regr parameters
  ini_list <- list(x0=matrix(1, mm*nn, 1))
  
  ## list of model matrices & vectors
  mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
  ##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
  
  ## list of control values
  con_list <- list(safe=TRUE, allow.degen=FALSE,
                   conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
  
  ## fit multivariate DLM
  mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
  
  ## get list of Kalman filter output
  kf_out = MARSSkfss(mod_res[[i]])
  
  ## ts of forecasted betas
  betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
  
  ## ts of SE of betas; nxT matrix
  betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
  
  ## ts of alphas
  alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
  ## ts of SE{alphas}
  alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
  
  ## ts of betas
  betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
  ## ts of SE{betas}
  betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
  
}
## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 2023 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2023 iterations. 
## Log-likelihood: -371.9797 
## AIC: 777.9595   AICc: 780.2954   
##  
##              Estimate
## R.diag       6.00e-01
## Q.alpha.var  5.16e-01
## Q.alpha.cov  5.16e-01
## Q.beta.var   6.42e-05
## Q.beta.cov   4.46e-05
## x0.X1       -1.94e-01
## x0.X2       -1.64e-01
## x0.X3       -5.67e-02
## x0.X4       -1.12e-01
## x0.X5       -2.22e-01
## x0.X6       -1.76e-01
## x0.X7        2.85e-01
## x0.X8       -7.16e-03
## x0.X9        1.94e-01
## x0.X10       2.58e-01
## x0.X11       1.92e-01
## x0.X12       1.46e-01
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 1943 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1943 iterations. 
## Log-likelihood: -378.4222 
## AIC: 790.8445   AICc: 793.2827   
##  
##             Estimate
## R.diag      6.29e-01
## Q.alpha.var 1.12e+00
## Q.alpha.cov 1.12e+00
## Q.beta.var  7.72e-05
## Q.beta.cov  5.71e-05
## x0.X1       1.42e-01
## x0.X2       1.23e-01
## x0.X3       1.12e-02
## x0.X4       6.52e-02
## x0.X5       1.45e-01
## x0.X6       9.14e-02
## x0.X7       3.97e-01
## x0.X8       5.01e-01
## x0.X9       4.26e-01
## x0.X10      4.81e-01
## x0.X11      3.22e-01
## x0.X12      2.90e-01
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 1688 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1688 iterations. 
## Log-likelihood: -140.0225 
## AIC: 302.0449   AICc: 304.1739   
##  
##              Estimate
## R.diag       2.15e-01
## Q.alpha.var  5.30e-01
## Q.alpha.cov  5.30e-01
## Q.beta.var   7.23e-05
## Q.beta.cov   5.88e-05
## x0.X1       -8.70e-01
## x0.X2       -8.96e-01
## x0.X3       -7.63e-01
## x0.X4        3.27e-01
## x0.X5        2.65e-01
## x0.X6        4.36e-01
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 1995 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1995 iterations. 
## Log-likelihood: -358.7276 
## AIC: 751.4551   AICc: 753.809   
##  
##             Estimate
## R.diag      4.92e-01
## Q.alpha.var 1.00e+00
## Q.alpha.cov 1.00e+00
## Q.beta.var  6.84e-05
## Q.beta.cov  5.03e-05
## x0.X1       2.22e-01
## x0.X2       1.16e-01
## x0.X3       1.29e-01
## x0.X4       2.23e-01
## x0.X5       1.86e-01
## x0.X6       9.38e-02
## x0.X7       5.22e-01
## x0.X8       1.97e-01
## x0.X9       3.74e-01
## x0.X10      6.10e-01
## x0.X11      4.37e-01
## x0.X12      6.92e-01
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop

3.7.1 Plots of results

betas <- betas[,names_mat$ID]

# total num of ts
NN <- dim(betas)[2]

# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
             rainbow(3, start=0.2, end=0.45, v=0.7),
             rainbow(6, start=0.8, end=0.95, v=0.85),
             rainbow(6, start=0.55, end=0.7))

# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7

# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))

# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)

# draw plots
for(i in 1:NN) {
    plot(t_index, betas[,i], type="n",
         ylim=yext, xaxt="n", yaxt="n",
         xlab="", ylab="", main=names_mat[i,"long"])
    rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
    abline(h=0, lty="dashed")
    box()
    lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
    lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
    lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
    #text(xin, yin, names.pop[i], adj=c(0,1))   
    if(i<=nr | nc==1) {
#       axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
#       axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
        axis(side=2, las=1, tick=TRUE)
        }
    if(i%%nr==0) {
        axis(side=1, at=t_index[(t_index %% 10)==0])
        }
    }
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)

3.8 v4: a diag & unequal; b diag & unequal

## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL

## number of regr params = level + slope
mm <- 2

## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
  
  ## print loop ID
  print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
  ## get popns
  assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
  ## number of assets
  nn <- dim(assets_mpg)[1]
  
  ## Process eqn
  ## B is Identity
  BB <- "identity"
  ## U is col vec of 0's
  UU <- "zero"
  ## Q will be block diagonal
  QQ <- matrix(list(0),mm*nn,mm*nn)
  ## upper L block for alpha
  aa <- 1:nn
  ## equal var-cov
  # QQ[aa,aa] <- "alpha.cov"
  # diag(QQ[aa,aa]) <- "alpha.var"
  ## diagonal & unequal
  diag(QQ[aa,aa]) <- paste0("alpha",seq(nn))
  ## lower R block for beta
  bb <- aa + nn
  ## diagonal & unequal
  diag(QQ[bb,bb]) <- paste0("beta",seq(nn))

  ## Observation eqn
  ## Z is array of regr variables (ie, 1's and market returns)
  ZZ <- array(NA, c(nn,mm*nn,TT))
  for(t in 1:TT) {
    ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
  }
  ## A is col vec of 0's
  AA <- "zero"
  ## assume same obs var & cov
  ##RR <- matrix(list(0),nn,nn)
  ##RR[aa,aa] <- "obs.cov"
  ##diag(RR[aa,aa]) <- "obs.var"
  RR <- "diagonal and equal"
  
  ## Initial states
  ##x0 <- rbind(matrix("pi.a",nn,1), matrix("pi.b",nn,1))
  ##V0 <- 100*diag(mm*nn)
  ##x0 <- rbind(matrix(0,nn,1), matrix(1,nn,1))
  ##V0 <- matrix(list(0),mm*nn,mm*nn)
  ##diag(V0[aa,aa]) <- "var.a"
  ##diag(V0[bb,bb]) <- "var.b"
  
  ## starting values for regr parameters
  ini_list <- list(x0=matrix(1, mm*nn, 1))
  
  ## list of model matrices & vectors
  mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
  ##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
  
  ## list of control values
  con_list <- list(safe=TRUE, allow.degen=FALSE,
                   conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
  
  ## fit multivariate DLM
  mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
  
  ## get list of Kalman filter output
  kf_out = MARSSkfss(mod_res[[i]])
  
  ## ts of forecasted betas
  betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
  
  ## ts of SE of betas; nxT matrix
  betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
  
  ## ts of alphas
  alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
  ## ts of SE{alphas}
  alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
  
  ## ts of betas
  betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
  ## ts of SE{betas}
  betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
  
}
## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 2398 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2398 iterations. 
## Log-likelihood: -423.8285 
## AIC: 897.657   AICc: 902.7751   
##  
##           Estimate
## R.diag    9.57e-01
## Q.alpha1  2.50e-05
## Q.alpha2  3.19e-01
## Q.alpha3  1.63e-01
## Q.alpha4  4.52e-05
## Q.alpha5  2.91e-05
## Q.alpha6  1.94e-02
## Q.beta1   1.81e-02
## Q.beta2   1.95e-04
## Q.beta3   4.10e-05
## Q.beta4   2.18e-05
## Q.beta5   3.96e-02
## Q.beta6   3.00e-05
## x0.X1    -2.03e-01
## x0.X2    -8.69e-01
## x0.X3     7.74e-01
## x0.X4    -3.23e-02
## x0.X5    -2.51e-01
## x0.X6    -2.97e-01
## x0.X7    -2.82e-01
## x0.X8    -1.24e-01
## x0.X9     1.22e-01
## x0.X10    1.80e-01
## x0.X11   -8.41e-01
## x0.X12    6.52e-02
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 1998 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1998 iterations. 
## Log-likelihood: -452.9958 
## AIC: 955.9916   AICc: 961.3414   
##  
##           Estimate
## R.diag    1.54e+00
## Q.alpha1  7.87e-05
## Q.alpha2  2.98e-01
## Q.alpha3  5.50e-05
## Q.alpha4  1.44e-04
## Q.alpha5  6.11e-05
## Q.alpha6  6.17e-05
## Q.beta1   7.04e-05
## Q.beta2   8.90e-05
## Q.beta3   2.55e-02
## Q.beta4   3.54e-03
## Q.beta5   9.93e-05
## Q.beta6   8.15e-05
## x0.X1    -1.17e-01
## x0.X2    -8.48e-01
## x0.X3    -3.39e-01
## x0.X4    -1.57e-01
## x0.X5    -1.55e-01
## x0.X6    -3.63e-01
## x0.X7     4.09e-01
## x0.X8     4.56e-01
## x0.X9    -8.80e-02
## x0.X10    3.41e-01
## x0.X11    1.81e-01
## x0.X12    3.54e-01
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 1257 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1257 iterations. 
## Log-likelihood: -177.6811 
## AIC: 381.3621   AICc: 384.3457   
##  
##           Estimate
## R.diag    2.24e-01
## Q.alpha1  6.78e-01
## Q.alpha2  2.49e-01
## Q.alpha3  3.72e-01
## Q.beta1   1.58e-04
## Q.beta2   5.87e-05
## Q.beta3   1.02e-04
## x0.X1    -3.95e-01
## x0.X2    -1.22e+00
## x0.X3    -5.32e-01
## x0.X4     2.65e-01
## x0.X5     2.16e-01
## x0.X6     4.79e-01
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 1405 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1405 iterations. 
## Log-likelihood: -458.2313 
## AIC: 966.4625   AICc: 971.6213   
##  
##           Estimate
## R.diag    5.73e-01
## Q.alpha1  7.26e-01
## Q.alpha2  4.67e-01
## Q.alpha3  2.71e-01
## Q.alpha4  3.54e-01
## Q.alpha5  8.82e-01
## Q.alpha6  1.66e+00
## Q.beta1   8.31e-05
## Q.beta2   6.73e-05
## Q.beta3   6.61e-05
## Q.beta4   5.74e-05
## Q.beta5   1.02e-04
## Q.beta6   2.26e-04
## x0.X1    -1.05e-01
## x0.X2    -4.65e-01
## x0.X3    -4.28e-01
## x0.X4     2.56e-01
## x0.X5     1.23e-01
## x0.X6     1.03e+00
## x0.X7     5.02e-01
## x0.X8     1.90e-01
## x0.X9     4.74e-01
## x0.X10    6.16e-01
## x0.X11    4.48e-01
## x0.X12    9.61e-01
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop

3.8.1 Plots of results

betas <- betas[,names_mat$ID]

# total num of ts
NN <- dim(betas)[2]

# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
             rainbow(3, start=0.2, end=0.45, v=0.7),
             rainbow(6, start=0.8, end=0.95, v=0.85),
             rainbow(6, start=0.55, end=0.7))

# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7

# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))

# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)

# draw plots
for(i in 1:NN) {
    plot(t_index, betas[,i], type="n",
         ylim=yext, xaxt="n", yaxt="n",
         xlab="", ylab="", main=names_mat[i,"long"])
    rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
    abline(h=0, lty="dashed")
    box()
    lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
    lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
    lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
    #text(xin, yin, names.pop[i], adj=c(0,1))   
    if(i<=nr | nc==1) {
#       axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
#       axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
        axis(side=2, las=1, tick=TRUE)
        }
    if(i%%nr==0) {
        axis(side=1, at=t_index[(t_index %% 10)==0])
        }
    }
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)